Official journal website: 
amphibian-reptile-conservation.org 


Amphibian & Reptile Conservation 
16(1) [General Section]: 94-105 (e305). 


Distribution and habitat suitability of two neighboring 
Lycian salamanders 


‘Omer Dilbe, ?Akin Kıraç, and **Eyup Başkale 


!Pamukkale University, Faculty of Science and Arts, Department of Biology, Denizli, TURKEY *Canakkale Onsekiz Mart University, Canakkale 
Technical Sciences Vocational College, Canakkale, TURKEY 


Abstract.—Lyciasalamandra fazilae and Lyciasalamandra flavimembris are two Endangered and endemic 
species which occur only in Mudla province of Turkey. In protecting an endemic or endangered species, the 
first step is to understand its potential and/or known distribution. Therefore, we used the Maximum Entropy 
modelling software (MaxEnt) to analyze the current potential distribution and most important habitat features 
associated with the localities of these two species. The variables with the highest contributions to the model 
were: Bedrock, Precipitation of Coldest Quarter, and Normalized Difference Vegetation Index for L. flavimembris; 
and Bedrock, Temperature Seasonality, Precipitation Seasonality, and Precipitation of Coldest Quarter for L. 


fazilae. We also identified two new localities for L. flavimembris using the habitat suitability model. 
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Introduction 


There are only seven species of Lycian salamanders in the 
world, six of which are found in Turkey. Among them, the 
Marmaris Salamander [Lyciasalamandra flavimembris 
(Mutz and Steinfertz 1995)] and the Gócek Salamander 
| Lyciasalamandra fazilae (Basoglu and Atatür, 1974)] 
are local endemic species distributed in the Muğla 
province of Turkey. Lyciasalamandra fazilae occurs 
in the eastern part of Mugla province (Fethiye, Gócek, 
Ortaca, and Köyceğiz districts), while L. flavimembris 
occurs in the western part of Muğla province (Milas, 
Ula, and Marmaris districts). Both species were formerly 
considered to be subspecies of Mertensiella luschani, 
with L. flavimembris even being con-subspecific with 
L. helverseni from the Greek Karpathos archipelago. 
However, previous studies have shown that they 
are morphologically and phylogenetically separate 
species (Oz et al. 2004; Veith et al. 2016, 2020; Veith 
and Steinfartz 2004), and their colorations are clearly 
distinguishable (Oz et al. 2004; Ozeti and Yilmaz 1994). 

Amphibians are highly susceptible to any changes 
in their habitat because of their highly permeable skin, 
and many species spend their lives in both terrestrial and 
freshwater habitats (Alford and Richards 1999; Barinaga 
1990; Duellman and Trueb 1994). The current information 
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on the ecology of Lycian salamanders broadly covers all 
species in this genus, and is therefore considered to be 
generally applicable to all of them (cf. Ozeti and Yilmaz 
1994). On this basis, the Lyciasalamandra species are 
terrestrial, inhabiting rocky limestone areas mostly 
in pine forests and maquis—sometimes near single- 
standing pines and olive trees, sometimes in deciduous 
forests dominated by oaks and junipers, and occasionally 
in accumulations of rocks or on slopes without vegetation 
(e.g., Baran and Atatur 1998; Başoğlu and Özeti 1973; 
Veith et al. 2001). The vertical distributions of these 
species are known to range from 25 to 1,400 m asl, 
where the mean annual rainfall may be less than 1,000 
mm (Veith et al. 2001; Yildiz and Akman 2015). 
Lyciasalamandra fazilae and L. flavimembris are 
listed as Endangered by the IUCN Red List of Threatened 
Species (http://www.iucnredlist.org; Accessed: 4 May 
2020) in view of their naturally restricted ranges and 
the continuing decline of their habitats. In protecting 
an endemic and/or Endangered species, the first step is 
to understand its potential and/or known distribution 
(Sousa-Silva et al. 2014). Intense research on the existing 
distribution of Lycian salamanders is time-consuming 
and expensive, but modelling their distributions could 
provide more accurate results with less time and effort 
(Hernandez et al. 2006). Species Distribution Modelling 
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Fig. 1. Study area and distributions of presence data for L. flavimembris and L. fazilae. 


(SDM) is a correlative approach in which habitat 
suitability, and therefore the distribution of a species, is 
estimated on the basis to environmental and geographical 
information (Elith and Graham 2009). The resulting 
models are called habitat suitability models, and they 
are considered to be important for the conservation of a 
species' habitat and the implementation of conservation 
action plans (Buckland and Elston 1993; Marzluff et 
al. 2002). They can be used to identify potential risks 
to a species and thus to prioritize habitat conservation, 
to optimize land management planning, and to allocate 
suitable habitats for potential translocation programs 
(Corsi et al. 1999; Ozkan and Berger 2014; Stoms et al. 
1992). 

The effective conservation of amphibian populations 
is typically limited by the lack of species-specific 
ecological knowledge. Therefore, this study was 
conducted to identify the environmental variables which 
limit the distribution of Marmaris Salamander and Gócek 
Salamander, and to determine their current and potential 
habitats. We believe that the models and maps obtained 
through the MaxEnt method will provide a base for the 
successful execution of species protection action plans. 


Materials and Methods 


Species data and study area. Between 2012 and 2020, 
field studies were carried out during the activity period 
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of the salamanders (October—April) within the province 
of Muğla, Turkey (Fig. 1). The study sites included 
four Specially Protected Areas (Gökova SPA, Datça- 
Bozburun SPA, Fethiye-Gócek SPA, and Köyceğiz 
Dalyan SPA), one National Park (Marmaris NP), and a 
Wildlife Development Area (Köyceğiz). The elevations 
of the sites ranged from 0 to 1,300 m asl. The climate 
is dominated by the Mediterranean climate. Urbanized 
areas, touristic areas, and natural areas without human 
intervention constitute important places in the study area 
which are mostly covered with maquis areas (shrublands), 
Red Pine (Pinus brutia) dominated coniferous forest, 
and agricultural fields. The field studies were carried out 
during both day and night. A total of 240 sample areas 
were examined, each with a size of 874 m x 874 m (Le., 
the resolution of the Worldclim [version 2.1] data used 
as described below). The altitudes and coordinates of 
each salamanders' presence point were recorded with a 
Garmin 62S GPS receiver using the WGS 84 coordinate 
system. 


Environmental data. The Aster Global Digital Elevation 
Model (GDEM), version 3, was obtained from Earthdata 
(http://earthdata.nasa.gov). Altitude, aspect, and slope 
were produced using GDEM (Zeiler 1999) in ArcMap 
10.2 software. The Topographic Position Index (TPI), 
Topographic Wetness Index (TWI), Landform Position 
Index (LPI), roughness index, hillshade index, ruggedness 
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index, solar radiation index, and solar illumination 
index (at 0600 h, 0800 h, 1000 h, 1200 h, 1400 h, 1600 
h, 1800 h, 2000 h, and total solar illumination) were 
created with the help of the “Topography tools" plugin 
included in ArcGIS 10.2 (Jenness 2006). The NDVI 
(Normalized Difference Vegetation Index) data produced 
by the MOD13Q1 module, which is one of the MODIS 
VI satellite data sources, was cut and used at the study 
area scale. NDVI values range from -1 to +1. Negative 
values represent water, snow, clouds, and non-plant 
areas; while positive values indicate the presence of 
vegetation. However, since negative values complicate 
the statistical analysis, the NDVI values were converted 
to the 0—10,000 range by using the formula: NDVI * 
10,000 (Celik and Gülersoy 2017). The bedrock map of 
the study area was obtained from the General Directorate 
of Mineral Research and Exploration (Maden Tetkik ve 
Arama Genel Müdürlüğü, http://yerbilimleri.mta.gov. 
tr/anasayfa.aspx). Different bedrock types (154) are 
shown in the form of polygons on the digital bedrock 
map obtained, which was used as a base map. These 
data were used as categorical data. Bioclimatic data 
representing the current climatic conditions of the study 
area were obtained from http://www.worldclim.org (Fick 
and Hijmans 2017). These data (Worldclim, Version 2.1) 
were obtained in the WGS 84 coordinate system with 
the highest resolution (30 arc-seconds, or 874 m x 874 
m), and in the ESRI Grid format. Nineteen bioclimatic 
variables (Biol—Biol9, Table 1) with this feature were 
cut on the scale of the study area with the help of ArcMap 
10.2. Temperature data (Biol, Bio2, Bio5—Bio11) values 
are shown multiplied by 100. 

For all of the digital base maps of the environmental 
variables in ASCII format, each cell was produced in the 
WGS 84 coordinate system (874 m x 874 m), and thus is 
of the same size as the sample areas. 


Statistical evaluation, habitat suitability model, and 
habitat suitability model map. Due to the small size of 
the study area, high correlation is expected between the 
bioclimatic data and other environmental variables. This 
may pose a problem during the analysis. To eliminate 
the multicollinearity problem, we applied Pearson 
Correlation Analysis, using a threshold of r?° < 0.8, for a 
total of 40 environmental variables. If a pair of variables 
was found to have a correlation coefficient greater than 
0.8, they were considered to represent related phenomena, 
and one of them was excluded from the analysis. 
Maximum Entropy (MaxEnt) (Phillips et al. 2006) is 
a popular habitat suitability modelling method, which 
provides more accurate results with less data in smaller 
areas compared with other methods (e.g.; DOMAIN, 
BIOCLIM, and GARP) (Hernandez et al. 2006; Phillips 
and Dudík 2008; Wisz et al. 2008). In addition, MaxEnt 
enables the joint processing of categorical and continuous 
data (Phillips and Dudík 2008), and it produces a habitat 
suitability map (Elith et al. 2011; Hernandez et al. 2006, 
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2008). MaxEnt is based on ENFA (Ecological Niche 
Factor Analysis; Hirzel et al. 2002) and examines the 
characteristics of the locations of the target species, and 
then estimates a suitability level for all areas based on 
the values taken by the factors which affect the known 
distribution of the species (Baldwin 2009). In this respect, 
the MaxEnt method was used to evaluate the potential 
distributions of Marmaris Salamander and Göcek 
Salamander using MaxEnt 3.4.1 software (Phillips et al. 
2006). MaxEnt calculates the maximum entropy to find 
the most likely geographical and ecological distribution of 
a target species. MaxEnt also examines the relationships 
between the asset data of the target species and 
environmental variables, and determines the ecological 
requirements of the target species. It then predicts the areas 
in which the target species will be more or less likely to 
appear based on the ecological requirements of the target 
species (Baldwin 2009). 

The environmental data, including presence data, in 
CSV format and environmental variables in ASCII format 
were analyzed with the help of MaxEnt 3.4.1 software. 
Species data were separated into 90% for training data 
and 1096 for test data using the software settings, and 
the analysis was adjusted to carry out ten repetitions. The 
replicated run type Crossvalidate was selected. Further 
settings were: maximum iterations — 500, convergence 
threshold = 0.00001, and default prevalence = 0.5. 
The Area Under the Receiver Operating Characteristic 
(ROC) Curves (AUC) was used to evaluate model 
performance. Finally, among the models with excellent 
model performance, the model with the lowest standard 
deviation between the training data AUC value and the 
test data AUC value was selected as “the best model,” 
and the species distribution maps of that model were 
visualized with ArcMap 10.2 software. 


Results 


For L. flavimembris and L. fazilae, 83 and 66 presence 
data points were obtained from the field studies, 
respectively, of which 68 and 54 were used for the final 
models, respectively. Most of the presence data obtained 
during the field studies were either known localities 
or points very close to known localities (Arslan et al. 
2018; Başkale et al. 2019; Göçmen et al. 2018; Oğuz 
et al. 2020; Polat and Baskale 2018; Veith et al. 2020). 
According to the results of the habitat suitability model, 
the training data set AUC value was 0.942 and the test 
data set AUC value was 0.941 + 0.056 (P « 0.001) for 
L. flavimembris (Fig 2a); while for L. fazilae the training 
data set AUC value was 0.954 and the test data set AUC 
value was 0.948 + 0.076 (P « 0.001) (Fig. 2b). These P 
values indicated that the model obtained was at the level 
of “perfect explanation" for the ecological requirements 
in the habitat preferences of both salamanders. 
According to the percentages of their contributions to 
the MaxEnt model, the important or highly contributing 
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(a) Sensitivity vs. 1-Specifity for L. flavimembris 
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(b) Sensitivity vs. 1-Specifity for L. fazilea 
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Fig. 2. The receiver operating characteristic (ROC) curves for L. flavimembris (a) and L. fazilae (b). 


variables which limit the geographical distribution ranges 
included three variables for L. flavimembris and four for 
L. fazilae (Fig. 3). The percentages of contribution to the 
model and occurrence intervals of these environmental 
variables are given in Table 1. The variables with the 
highest contributions to the model were: Bedrock, 
Precipitation of Coldest Quarter, and Normalized 
Difference Vegetation Index for L. flavimembris (Table 
1 and Fig. 4a); and Bedrock, Precipitation of Coldest 
Quarter, Temperature Seasonality, and Precipitation 
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Seasonality for L. fazilae (Table 1 and Fig. 4b). Combined, 
these variables explained 86.3% and 99.1% of the 
variation in the two species distributions, respectively. 
The habitat suitability models showed the potential 
distributions of the two species, and the predicted models 
confirmed the mostly known geographical ranges of both 
of them (Fig. 5). The area of high predicted probability of 
occurrence for L. flavimembris was concentrated around 
the Kótekli, Ula, Milas, and Marmaris districts (Fig. 5a). 
In particular, the southwestern part of Marmaris district is 
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Fig. 3. Results of the Jackknife test for evaluating the relative importance of environmental variables for L. flavimembris (a) and L. 
fazilae (b). See Table 1 for definitions of the environmental variables. 


the most intensely occupied area for L. flavimembris. In 
relation to the habitat suitability model of L. flavimembris, 
the field studies revealed two new localities: Kizilkóy 
(36?41'N, 28°06’ E; 204 m asl) in the Selimiye district, 
and İçmeler (36°46’N, 28?12'E; 142 m asl) in the 
Marmaris district. For L. fazilae, the habitat suitability 
model indicated a high probability of occurrence mostly in 
known habitats, such as Gókgeovacik (Fethiye), Üzümlü 
(Fethiye), Dalyan (Ortaca), Kapikargin (Dalaman), and 
Sultaniye (Köyceğiz) (Fig. 5b). 


Discussion 


According to the MaxEnt results, the average 
contributions (in percentage) of the key environmental 
variables to the model were determined as: Bedrock 
(54.2%), Precipitation of Coldest Quarter (41.7%), and 
NDVI (4.1%) for L. flavimembris, and Bedrock (62.3%), 
Precipitation of Coldest Quarter (25.4%), Temperature 
Seasonality (8.5%), and Precipitation Seasonality (3.8%) 
for L. fazilae. 

The most important factor limiting the distributions 
of both L. flavimembris and L. fazilae is bedrock type 
rather than any of the climatic conditions. Species that 
prefer specific bedrock types need corridors made up of 
suitable bedrock to expand their distributions (Sinervo 
et al. 2017). It is known that salamanders which live in 
suitable bedrock often hide in the cracks, cavities, and 
underground of this bedrock under unfavorable climatic 
conditions (Baran and Atatur 1998). These cracks 
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and holes maintain proper moisture and temperature 
conditions. The MaxEnt outputs of environmental 
variables showed that L. flavimembris prefers 10 of the 
154 bedrock types in the region, while L. fazilae prefers 
nine bedrock types (Fig. 6). While previous studies 
revealed only limestone (Göçmen and Karis 2017; Veith 
et al. 2001), this study shows that L. flavimembris and 
L. fazilae can be found under different types of stones 
but their habitats mostly include limestone and cherty 
limestone. 

Amphibians have a high climatic sensitivity. due 
to their ectothermic physiology and their constant 
need for moisture (Wells 2007). Previous studies have 
emphasized that humid areas, areas with a dense green 
cover, an average annual rainfall of 800—1,500 mm, and 
rocks with moist ground crevices are suitable habitats 
for Lycian salamanders (Baran and Atatur 1998; Veith et 
al. 2001). Ródder et al. (2011) investigated the climatic 
niche similarities between the Lycian salamander species 
using 19 bioclimatic data sets. That study found that 
Lycian salamanders (except for L. helverseni) preferred 
similar climatic conditions, and the mean Temperature of 
Coldest Quarter (variable Bio11) ranged from 6-12.5 °C 
and Precipitation of Coldest Quarter (B1019) ranged from 
350—620 mm. The species-specific studies have shown 
that Pinus brutia, Mediterranean maquis, green mosses, 
and limestones are indicators for L. flavimembris habitat 
(Göçmen and Karis 2017), and the air temperature interval 
of the active season of L. flavimembris ranged from 5-21 
°C, while monthly average precipitation ranged from 
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Fig. 5. MaxEnt habitat suitability maps for L. flavimembris (a) and L. fazilae (b). 
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Dolomite [Dol], Limestone [LS], Melange [Mel], Peridotite [Per], Spilite-Basalt-Tuff [SBT], Sandstone-Mudstone [SS-MS], 
Sandstone-Mudstone-Limestone [SS-MS-LS], Volcanite-Sedimentary Rock [V-SR], and all unknown rock types [Unknown]. 


57—335 mm (Baskale et al. 2019). On the other hand, 
Polat and Baskale (2018) stated that the greatest number 
of individuals of L. fazilae was observed at temperatures 
between 2 and 18 °C (mean 12.99 + 0.403 °C), and that 
the active period started with the first autumn rains and 
a sharp decrease in air temperature (« 20 ?C), and ended 
with higher air temperatures (22 °C and above). 

Climatic conditions may also limit the distributions 
of both species, resulting in narrow distribution areas. 
Our habitat suitability models show that L. flavimembris 
and L. fazilae both have specific demands with respect 
to precipitation and temperature. Specifically, the 
Precipitation of Coldest Quarter (Bio19) is 500-650 mm 
for L. flavimembris, while for L. fazilae the Precipitation 
of Coldest Quarter (Biol9; 600 mm), Precipitation 
Seasonality (Biol5; 100 mm), and Temperature 
Seasonality (Bio4; 5—6 ?C) were found to be important 
predictors of its distribution. These results show that the 
current climatic conditions are sufficient for L. fazilae 
and L. flavimembris to survive. This supports the MaxEnt 
ClogLog values for L. fazilae and L. flavimembris given 
in Veith et al. (2020), which showed the prevalence of 
unsuitable current climatic conditions for the survival of 
many of the Lycian salamanders other than L. fazilae and 
L. flavimembris. 

In our models, vegetation is another of the 
environmental factors that determine the distributions 
of the two salamander species. Lyciasalamandra 
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flavimembris was detected in areas with an NDVI of 
10,000, indicating green areas with high canopy cover. 
For L. fazilae, the interval of the NDVI value was 
wider (1,000-10,000), hence its distribution area is 
characterized by more heterogeneous vegetation, such 
as pine forests, Mediterranean marquis, and olive tree 
fields. Our habitat compatibility model obtained with the 
MaxEnt method is compatible with the known biology 
of Lycian salamanders (Baran and Atatur 1998; Ozeti 
and Yilmaz 1994; Veith et al. 2001). Another consistency 
in our results is that the locations with the highest 
population densities and abundances of L. flavimembris 
and L. fazilae shown in Polat and Baskale (2018) and 
Baskale et al. (2019) are the same as the localities with 
high suitability values in our habitat suitability map. 

Our habitat suitability maps mostly reflect the known 
localities of both species, but it is important to consider 
some differences between the predicted model and 
the known habitats. For L. flavimembris, the habitat 
suitability map shows inhabitable areas to the west. 
Although the Yalikavak and Mazi Mountain (Milas) 
populations (Oguz et al. 2020) are located in this area, 
the potential distribution is extended even to the Bodrum 
district. This suggests that there are either important 
barriers to the species’ dispersion, or it has simply not 
yet been recorded from these areas. Moreover, the 
model predicted suitable habitats for L. fazilae within 
the distribution area of L. flavimembris (see also Veith 
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et al. 2020). This situation arises from the fact that both 
species prefer similar environmental variables such as 
Precipitation of Coldest Quarter and Bedrock. On the 
other hand, Veith et al. (2020) showed a strong degree of 
isolation among Lyciasalamandra populations, including 
phyloclades of L. fazilae, and two subspecies of L. fazilae 
are recognized: L. f. fazilae and L. f. ulfetae (Göçmen 
et al. 2018). However, Veith et al. (2020) claimed that 
the L. fazilae phyloclade diversity is higher than that 
reflected by current taxonomy, with five phyloclades 
forming three well-supported phylogenetic clusters: 
(faz-I + faz-II), faz-III, and (faz-IV + faz-V). The vertical 
extension of the Taurus Mountains between the Gócek 
and Dalaman districts constitutes the first (faz-I + faz-II) 
and the second (faz-III) phylogenetic clusters. However, 
the third cluster (Ülemez population and Sultaniye 
population) is geographically isolated by the Köyceğiz 
Lake and Dalyan Canal in the east, and the Ulemez 
Mountain and the extensions of Taurus Mountains in the 
west and northwest (see Figs. 1 and 5b). 

In conclusion, potential distribution maps of L. 
flavimembris and L. fazilae were created based on 
bioclimatic data and some environmental variables. These 
maps indicated that the current climatic conditions of 
the regions where both species live are suitable for the 
survival of the species. In addition, some populations of L. 
flavimembris (1.e., Yaylasógüt and Arıcılar) and L. fazilae 
(i.e., Üzümlü and Gékceovacik) were located far from 
the Mediterranean coast, indicating that these species can 
tolerate more diverse climatic conditions. This study is an 
important step for the conservation of endangered species 
within and outside existing protected areas, and may help 
alleviate the population decline of both species. 
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